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We study the equation of state at finite temperature and density in two-flavor QCD with the RG- 
improved gluon action and the clover-improved Wilson quark action on a 16"^ x 4 lattice. Along the 
lines of constant physics at mps/my = 0.65 and 0.80, we compute the second and forth derivatives of 
the grand canonical partition function with respect to the quark chemical potential Hq — (/i„ + ^td)/2 
and the isospin chemical potential fij = {fiu — l-i'd)/2 at vanishing chemical potentials, and study 
the behaviors of thermodynamic quantities at finite fiq using these derivatives for the case fii — 0. 
In particular, we study density fluctuations at nonezero temperature and density by calculating 
the quark number and isospin susceptibilities and their derivatives with respect to fiq. To suppress 
statistical fluctuations, we also examine new techniques applicable at low densities. We flnd a 
large enhancement in the fluctuation of quark number when the density increased near the pseudo- 
critical temperature, suggesting a critical point at flnite fiq terminating the flrst order transition 
line between hadronic and quark gluon plasma phases. This result agrees with the previous results 
using staggered-type quark actions qualitatively. Furthermore, we study heavy-quark free energies 
and Debye screening masses at finite density by measuring the first and second derivatives of these 
quantities for various color channels of heavy quark-quark and quark-anti-quark pairs. The results 
suggest that, to the leading order of fiq, the interaction between two quarks becomes stronger at 
flnite densities, while that between quark and anti-quark becomes weaker. 

PACS numbers: 11. 15. Ha, 12.38.Gc, 12.38.Mh 

I. INTRODUCTION 

Heavy-ion collision experiments are taking place at BNL aiming at the experimental studies of a new state of 
matter, the quark-gluon plasma In order to extract unambiguous signals for the QCD phase transition from the 
heavy-ion collision experiments, quantitative calculations directly from the first principles of QCD are indispensable. 
At present, the lattice QCD simulation is the only systematic method to do so. Various computational techniques have 
been developed to study the nature of quark matter at finite temperature (T) and at small chemical potentials fj,u and 
[IjQ- From intensive studies for the isosymmetric case Hu = l^d = f^q, it turned out that accurate zero-temperature 
simulations are important to set the scale to achieve high precision results at finite T and fiq. 

Most of the lattice QCD studies at finite so far have been performed using staggered-type quark actions with 
the fourth-root trick for the quark determinant. However, the fourth-root trick makes the theory non-local and thus 
the universality arguments fragile. It should be kept in mind that the staggered-type quarks for two-favor QCD does 
not show the scaling properties at finite T expected from the three-dimensional 0(4) spin model [1, @. This may 
suggest large lattice artifacts to the results of staggered-type quarks near the transition point. Moreover, problems in 
the staggered quark formulation at finite density are pointed out in ^ . Since the theoretical base for the fourth-root 
trick is not clear, it is indispensable to carry out simulations adopting different lattice quark actions to control and 
estimate systematic errors due to the lattice discretization. 

Several years ago, the CP-PACS Collaboration has studied finite-temperature QCD using the clover-improved 
Wilson quark action coupled with the RG-improved Iwasaki action for gluons 0, 13 . With two flavors of dynamical 
quarks, the phase structure, the transition temperature and the equation of state have been investigated. In contrast 
to the case of the staggered-type quarks, both the standard Wilson quark action [91] and the clover-improved Wilson 
quark action Q reproduce the expected universality around the critical point of the chiral phase transition: the 
subtracted chiral condensate shows the scaling behavior with the critical exponents and scaling function of the three- 
dimensional 0(4) spin model. Moreover, extensive calculations of rnajor physical quantities such as the light hadron 
masses have been carried out at T = using the same action [lOi. ill|. Therefore, it is worth revisiting this action 
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armed with recent techniques for finite fiq. 

In the (T, jig) plane, phenomenological studies suggest the existence of a critical point at which the first order 
phase transition line separating the hadronic phase and the quark-gluon-plasma phase terminates [l2l - [l^ . Because 
the critical point has second order characteristics, the fluctuation of the net quark number will diverge as we approach 
to the critical point in the (T, /ig) plane, while the fluctuation in the isospin number will remain finite [isl [la|. Such 
hadronic fluctuations may be experimentally examined in heavy-ion collisions by an event-by-event analysis. The 
Bielefeld- Swansea Collaboration reported lattice results for the quark number susceptibility (the second derivative of 
the thermodynamic grand canonical potential lu/T"^ = —{VT^)~^\nZ, which is proportional to the pressure of the 
system) by the Taylor expansion method using a p4-improved staggered quark action p7l - [l9| : From a calculation of 
the Taylor expansion coefficients of w/T'* up to 0[{fiq/T)^], they found that the quark number fluctuation increases 
rapidly as fig increases in the region near the transition temperature. This suggests indirectly the existence of the 
nearby critical point in the (T,fiq) plane. Moreover, 2-I-I flavor simulations in staggered quarks with almost physical 
quark masses have recently been performed and the same behaviors in the fluctuations have been found at finite 
density [13, [2l|. Therefore, it is important to confirm the result using the Wilson- type quarks. 

In this paper, we study thermodynamic properties of QCD at finite temperature and density with two flavors of 
clover-improved Wilson quarks coupled with the RG-improved Iwasaki gluons. The simulations are performed along 
lines of constant physics corresponding to the pion and rho meson mass ratio, mps/my = 0.65 and 0.80 at T = 0. 
We calculate the Taylor coefhcients for the pressure in terms of Hq/T up to the fourth order, and study the quark 
number and isospin susceptibilities at finite /iq. Since the odd derivatives vanish at /i^ — 0, the fourth derivative 
is the leading contribution to the /i^-dependence of susceptibilities. We find that Wilson-type quarks require much 
more statistics than staggered-type quarks to obtain the susceptibilities with a comparable quality. To overcome this 
problem, we introduce a couple of tricks in the evaluation of the Taylor expansion coefficients. Furthermore, we adopt 
a hybrid method of Taylor expansion and spectral reweighting in which w/T^ for the reweighting is approximated 
by a truncated Taylor expansion p^ . [2^ . Since the applicable range of the reweighting method is narrow due to the 
sign problem, we introduce the Gaussian method proposed in [2j|. Using these techniques, we compute the quark 
number density and the susceptibility in a relatively wide range of fJ-q/T, and compare the results with those with 
staggered-type quarks. 

We also extend our previous study of heavy-quark free energies in various color channels at fiq — (23 | to finite 
lj,q. At T > Tpc, where Tpc is the pseudo-critical temperature, we calculate the Taylor expansion coefficients for 
the heavy-quark free energies between a static quark (Q) and an antiquark (Q) and those between Q and Q, for all 
color channels up to the second order in ^iq/T . By comparing the expansion coefficients of the free energies, we find 
that the inter-quark interaction between Q and Q becomes weaker, whereas that between Q and Q becomes stronger 
as pLq increases. The expansion coefficients of the effective running coupling acff(r, /^g) and the Debye screening 
mass mD(T, iiq) are also extracted by fitting the numerical results with a screened Coulomb form; we find that the 
heavy-quark free energies are well reproduced by the channel dependent Casimir factor and the channel independent 
aes{T,fJ,q) and mD{T,^q) at T ^ 2rpc- Magnitude of the second order coefficient of mu{T,^q) does not agree with 
that of the leading-order calculation in the thermal perturbation theory. 

In Sec. [ni we summarize our lattice action and simulation parameters, and determine the pseudo-critical temper- 
ature. In Sec. mil we calculate the Taylor expansion coefficients of the therniodinamic grand canonical potenital in 
terms of the quark chemical potentials and /i^ and evaluate them for the isosymmetric case Hu = jJ-d = at 
/j,g = up to 0(/j,g). In Sec. llVi we adopt the hybrid method combined with the Gaussian method, to improve the 
calculation. The static quark free energies and the Debye screening mass are discussed in Sec. El Conclusions and 
discussions are given in Sec. IVII We summarize properties of the pressure and the quark number susceptibility in 
the free gas limit in Appendix [X] Appendix [B] is devoted to a description of detailed derivations of formulae for the 
Gaussian method. Results of the fits of heavy-quark free energies are summarized in Appendix [C] 

II. PHASE STRUCTURE AND LINES OF CONSTANT PHYSICS AT ^i, = 

A. Lattice action 

First, we summarize our simulation details. We adopt the same lattice actions as in our previous study at /Ug = 
[2^ . We use the RG-improved Iwasaki gauge action [2^ and the Nf ^ 2 clover- improved Wilson quark action [2^ 
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FIG. 1: Lines of constant physics (LCP) determined by mps/mv at T = (solid lines) for mps/my ~ 0.65, 0.70, 0.75, 0.80, 
0.85, 0.90 and 0.95. Kc is the chiral limit, i.e. mps/mv ~ 0. Dashed lines represent lines of constant T/Tpc on Nt = 4: lattices, 
where Tpc is the pseudo-critical temperature corresponding to Kt{Nt — 4) shown by the thick dashed line. 



defined by 



S — Sg -\' Sqj 

f=u,d x.y 



(1) 

(2) 
(3) 



where /3 — ci — —0.331, cq = 1 — 8ci and 



1=1 

-K{e''{l - lA)Ua,A5^+i,y + e-^{l + lAWl^S^ y^-^} - S^yCswK ^ a^,F^,. 



(4) 



Here K is the hopping parameter, /i = /x^a is the quark chemical potential in lattice unit and F^^, is the lattice 
field strength, F^^, = {ff^^, — fl^)/{8i), with the standard clover-shaped combination of gauge links. For the 
clover coefficient csw, we adopt a mean field value using Vt^^^^ calculated in the one- loop perturbation theory [25j : 
csw — (W^^^)^^^^ = (1 — 0.8412/9"^)"^/'*. We denote the spatial and temporal lattice size as Ns and Nt respectively. 
At Uq = 0, the phase diagram of this action in the (/3, K) plane has been obtained by the CP-PACS Collaboration 
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For phenomenological applications, we need to investigate the temperature dependence of thermodynamic observ- 
ables in a given physical system. On the lattice, "a given physical system" corresponds to a given set of values of 
dimension-less ratios of physical observables at T = and fiq = 0. Assuming the scaling, this forms a line in the 
coupling parameter space, called the line of constant physics (LCP), along which the lattice scale (lattice spacing a) 
is varied for a given physical system. On a finite-temperature lattice with fixed Nt, the temperature, T = 1/Nta, 
is varied along a LCP according to the variation of a. In this study, we determine LCP by mpg/my (the ratio of 
pseudo-scalar and vector meson masses at T = and /ig = 0). The bold solid line denoted as Kc in Fig. [T] represents 
the chiral limit, i.e. mps/wv = 0- Above the Kc fine, the parity-flavor symmetry is spontaneously broken 27]. 
The region below Kc corresponds to the two-flavor QCD with finite quark mass. We perform simulations in this 
region. The lines of constant mps/mv are investigated in Refs. [1, [2^, which is shown as thin solid lines in Fig. [1] 
corresponding to mps/my = 0.65, 0.70, 0.75, 0.80, 0.85, 0.90 and 0.95. 

The temperature T is estimated by the zero-temperature vector meson mass myalP, K) using 



T 1 

rrw^'^'^^ ~ Nt X mya{^, K) ' 



(5) 
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TABLE I: Simulation parameters for nips/niY = 0.65 (left) and mps/my = 0.80 (right) on a 16'^ x 4 lattice. 
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1.07(4) 


5000 


1 


85 





143502 
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0.80(4) 


6000 


1.70 





142871 
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FIG. 2: if-dependence of Polyakov loop for Nt = 4 (left) and 6 (right). Data at /3 = 1.7 and 1.8 for iVt = 4 and 1.9 and 1.95 
for Nt = 6 are renewed from Refs. 0, 



The lines of constant T/Tpc is determined by the ratio of T/my to Tpc/my where Tpc/my is obtained by T /my at 
Kt on the same LCP. We use an interpolation function, Tpc/my = A{\ + B{mps/my)'^)/{l + C{mps/my)^) with 
A = 0.2253(71), B = -0.933(17) and C = -0.820(39), obtained in Ref. @ to evaluate Tpc/my for each mps/my. 
The bold dashed line denoted as Kt{Nt — 4) in Fig. [T] represents the pseudo-critical line T/Tpc ~ 1. The thin dashed 
lines represent the results for T/Tp^ = 0.8, 1.2, 1.4, 1.6, 1.8, 2.0 at Nt = 4. 

We perform finite temperature simulations on a lattice with a temporal extent Nt = 4: and a spatial extent Ns — 16 
along the LCP's at TOps/my — 0.65 and 0.80. The standard hybrid Monte Carlo algorithm is employed to generate 
full QCD configurations with two flavors of dynamical quarks. The length of one trajectory is unity and the step 
size of the molecular dynamics is tuned to achieve an acceptance rate greater than 70%. Runs are carried out in the 
range /3 = 1.50-2.40 at thirteen values of T/Tpc ^ 0.82-4.0 for mps/my = 0.65 and twelve values of T/Tpc 0.76-3.0 
for mps/my = 0.80. Our simulation parameters and the corresponding temperatures are summarized in Table ID 
Because the determination of the pseudo critical line is more difficult than the calculation of T/my, the dominant 
source for the error of T/Tpc in Table |T] is the overall factor Tpc/my. The number of trajectories for each run after 
thermalization is 5000-6000. We measure physical quantities at every 10 trajectories. The study of heavy quark free 



energies a.t — using the same configurations have been already published in Ref. [24 1 



B. Critical temperature 

We update the analysis of the pseudo critical temperature done in Refs. 0, performing additional simulations 
at /3 = 6/g^ = 1.7 and 1.8 on an x iV* = 16 x 4 lattice and at 1.9 and 1.95 on x Nt ^ 16^ x 6. The number 
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TABLE II: Finite temperature transition/crossover point Kt for Nt = 'i and 6. Results for mps(r = 0)/mv(T = 0), mpsa(r : 



0), 



a{T > 0), Tpc/mY{T = 0), Tpc/y/^, Tpcro, and mpsro, are interpolated to the Kt line, 



mps/mv mpsa 



_AWI — 
rrig a 



Tpc/my Tpcj^fa Tpcrp 



mpsro 



Nt 



1.700 0.15014(33) 

1.800 0.14425(16) 

1.850 0.14019(18) 

1.900 0.13621(15) 

1.925 0.13417(23) 

1.950 0.13040(97) 

2.000 0.12371(73) 

2.100 0.10921(43) 



0.151987(22) 
0.147678(15) 
0.145526(58) 
0.143737(48) 

0.142072(14) 
0.140811(55) 
0.139020(21) 



0.509(35) 
0.7070(79) 
0.7905(60) 
0.8525(39) 

0.9051(64) 
0.9450(36) 
0.9790(13) 



0.579(51) 
0.849(18) 
1.031(15) 
1.183(11) 

1.440(66) 
1.689(39) 
2.196(18) 



0.1107(77) 

0.1864(72) 

0.2464(49) 

0.2725(67) 

0.363(25) 

0.500(18) 



0.2197(47) 

0.2083(21) 0.4204(29) 0.4716(42) 1.601(37) 

0.1917(20) 0.4359(60) 0.484(11) 1.994(55) 

0.1801(12) 0.4382(70) 0.484(16) 2.290(79) 

0.1572(62) 
0.1398(29) 
0.1114(9) 



iVt = 6 



1.950 0.14090(13) 0.142072(14) 0.591(21) 0.448(24) 0.0451(51) 

2.000 0.13861(21) 0.140811(55) 0.725(16) 0.580(27) 0.080(10) 

2.100 0.13365(40) 0.139020(21) 0.8635(78) 0.821(34) 0.175(13) 

2.200 0.12539(25) 0.137658(53) 0.9481(19) 1.240(16) 0.3607(67) 

2.300 0.11963(15) 0.136513(85) 0.9724(12) 1.454(8) 0.4813(39) 



0.2202(44) 0.4336(40) 0.4973(58) 1.336(73) 
0.2086(53) 0.4639(77) 0.530(13) 1.842(98) 
0.1753(58) 0.491(12) 0.570(13) 2.81(13) 
0.1275(15) 

0.1114(6) 



of trajectories for each new run is 1050-4200 after thermalization. We add the new data to the data in Refs. 0, 
and determine the pseudo-critical hopping parameters Kt defined from the peak of the Polyakov loop susceptibihty 
on 16'^ X 4 and 16"^ x 6 lattices, as a function of /3. Figures [5] and [3] are the results of the Polyakov loop {L) and 
Polyakov loop susceptibility xli respectively. We find a pronounced peak in the Polyakov loop susceptibility except 
for (3 — 1.90 at Nt = 6. The peak position of the susceptibility (Kt) is determined by fitting three or four data near 
the peak with the Gaussian form. The results are summarized in Table together with values of some quantities at 
Kt to set a physical scale. 

We use the data of the pseudo-scalar and vector meson masses at T = 0, wps and my, summarized in the Table 
IV of Ref. [1], and interpolate them following the method discussed in Refs. 01 and We also calculate the current 
quark mass defined through an axial vector Ward-Takahashi identity, V^yly^ = 2m_^^P + 0{a), where P is the 
pseudo-scalar density and Af^ the /z-th component of the local axial vector current [H, [2^ . Because the T-dependence 
in is small, we use the data of obtained in finite temperature simulations at iVj = 4 and 6 0, Q ■ In Table 

im on the Kt line are obtained using a cubic spline interpolation for each /3. A straight line interpolation leads 

to almost identical results within statistical errors. The values of the string tension a and the Sommer scale Tq [30j 
are estimated by interpolating or extrapolating the data at /3 = 1.80, 1.95, 2.10 and 2.20 TT] in the {j3,\/K — ^/kJ) 
parameter plane. 

The results of the pseudo-critical temperature are also shown in Table |TT1 We plot Tpc/my as a function of 
(mps/mv)^ in Fig. [H and find that the results of TVt = 4 and 6 agree with each other. Note that Tpc/my vanishes in 
the heavy quark limit mpg/mv — 1. Figure [J] suggests Tpc/my ~ 0.22 (Tpc ~ 170MeV) in the chiral limit. 

We denote the critical temperature in the chiral limit as Tc- As discussed in 0, the subtracted chiral condensate 
[29I satisfies the scaling behavior with the critical exponents and scaling function of the 3-dimensional 0(4) spin 
model. For the reduced temperature t and external magnetic field h, we adopt t ~ /3 — /Set and h ~ m^, where /Set 
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TABLE III: The critical point (/3ct) and critical temperature (Tc) in the chiral limit obtained by various fitting procedures. The 
fit range for /3 is written in "/3 range". Tc in a physical unit is estimated from the vector meson mass my = nip — 770MeV. 
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FIG. 4: Tpc/mv vs. mps/mv for Nt = 4: (circle) and 6 (triangle). The lightest two points for A'^t = 4 and the lightest one point 
for Nt = 6 are updated from Ref. Q. 



is the critical transition point in the chiral limit. For a precise determination of Tc, we need to deduce /3ct from 
the data. In this study, we perform critical scaling fits assuming that the pseudo-critical temperature tpc from the 
Polyakov loop susceptibility, as well as that from the chiral condensate, follows the scaling law tpc ~ h^^'^ with the 
0(4) critical exponent 1/y = 1/{PS) ~ 0.537(7). In practice, we fit the data of f3pc{K), i.e. the inverse function of 
KtiP) in Table m by 



I3pc = Pet + Ah^'y 



(6) 



with two free parameters, Pet and A. 

For the quark mass niq ^ h in the scaling fits, we test three variants. The first is m^a ~ l/X — 1/-K^c, where 
is the chiral point at which the pion mass vanishes at T = for each /3. The second is niqa ^ (mpga)^. The third 
is m^^a, i.e. the quark mass defined by the axial vector Ward-Takahashi identity. We plot (ipc as a function of 
\/K — 1/Kc (left), (mpsa)^ (center) and rriq^^a (right) in Fig. [5l The results of (ict and Tc are summarized in Table 
mil where Tc in MeV is calculated by Tc = \/[Nta{(ict)] with a from the vector meson mass mY{T = 0) = rup = 770 
MeV at Pet on Kc- We test two fit ranges of P for each extrapolation, which is denoted in Table Hill as "/3 range". We 
note that these 0(4) fits reproduce the data of Ppc much better than a naive linear fit Ppc — Pet + Ah. A tentative 
conclusion is that the critical temperature in the chiral limit is in the range 171-180 MeV for Nt — A and 160-184 
MeV for Nt — 6. There is still a large uncertainty from the choice of the fit ansatz and the fit range. To remove this, 
further simulations at lighter quark masses are necessary. 

For a comparison with other groups, we estimate Te in units of the Sommer scale rg [30] at Pet in the chiral limit for 
Nt = 6. Using the data of rp/a in the chiral limit at = 1.80, 1.95 and 2.10 [lH, we interpolate a/ro by a quadratic 
function and calculate TcTq — {Nta/ro)~^ . The estimates are about Tcrp « 0.40, as listed in Table Hill These values 
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FIG. 5: The pseudo-critical point /3pc as a function of [mqaY^'" with niqa ~ 1/K-1/Kc{\eit), (m|sa)^ (center) and m^^'a(right) 
for A''t = 4 (circle) and Nt = 6 (square). We fit the data in two fit ranges. The solid and dashed lines are the fit results with 
the long and short fit ranges, respectively. 



are close to TcVq = 0.402(29) obtained by the MILC Collaboration using the asqtad quark action in 2+1 flavor QCD 
[U ^ On the other hand, the RBC-Bielefeld Cohaboration obtained T^ro = 0.444(6)^3^ using a 2+1 flavor p4fat3 
improved staggered quark action 3,^. From a simulation of 2 flavor QCD using a clover improved Wilson action and 
the standard one-plaquette gauge action, the DIK Collaboration obtained T^ro = 0.438(6)j^7^ at the physical pion 
mass point, and the value in the chiral limit is 2% smaller than this value [33|. Our result is somewhat smaller than 
these values. Finally, the Budapest- Wuppertal group used a stout-link improved staggered fermion action and fixed 
the scale by the pion decay constant /jr. They found that Tc determined by the chiral susceptibility is Tc — 151(3)(3) 
MeV and that by the renormalized Polyakov loop is Tc = 176(3) (4) MeV in the continuum limit at the physical point 
[34j . Our result is close to their result defined by the Polyakov loop. For further discussions, see Refs. (35l - [37| . 



III. EQUATION OF STATE AT FINITE DENSITIES BY THE TAYLOR EXPANSION METHOD 

The main difficulty in a study of QCD at finite density is that the Boltzmann weight is complex for nonzero fiq. 
The quark matrix at zero density have the 75 Hermiticity = 75M75 which guarantees that the quark determinant 
is real. However, at /iq ^ 0, we have only 

Mt(/i,)=75Af(-^,)75, (7) 

from Eq. Therefore, the quark determinant is complex for 7^ 0. 

Because configurations cannot be generated with a complex probability, the conventional Monte Carlo method 
is not applicable at fig ^ 0. At present, there are three methods to study finite density QCD, all of which are 
applicable for small lig regions. The simplest is the method based on a Taylor expansion in terms of fJ^q/T around 
Mg = [Hi m, ■ Because the simulations at /ig = is free from the complex weight problem, the expansion 
coefficients, i.e. derivatives of physical quantities with respect to iig/T, can be evaluated by a conventional Monte 
Carlo simulation. The second approach is the reweighting method |40l - t43| . Performing simulations at Hq = 0, 
expectation values at finite are computed adopting a corrected Boltzmann weight. For the correction, quark 
determinant at finite /i^ is estimated numerically. Because fiuctuations in the complex phase of the determinant are 
large at large ^q and/or large lattice volume, a reliable calculation of expectation value becomes gradually difficult 
off the small /iq and small lattice volume region due to the sign problenr [4^ l45j . The third approach is the analytic 
continuation from simulations with imaginary chemical potentials j46l |47||. Since the equation ([7]) is generalized to 
M'l'(yLtg) — ■'j^M (— Aiq)75 for complex fiq, the Boltzmann weight is real and simulations are possible when the chemical 
potential is purely imaginary. Using results by the imaginary chemical potential simulations, information at a real 
chemical potential can be obtained by an analytic continuation. The analytic continuation is usually based on a 
Taylor expansion in terms of /i^ around fiq = for the study in the low density region, and improvements of the 
analytic continuation have been also discussed in [isl - lsoj to obtain reliable results in a wide range of real ^q. 



^ Originally, Tc is given in units of ri in Ref. [3l| . The scale of has been converted to rg using rg/ri = 1.4795 [3^ . 
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In this section, we adopt the Taylor expansion method to study the effects of /iq in the equation of state. Most of 
thermodynamic quantities, such as energy density, quark number, order parameters and various susceptibihties, are 
given by derivatives of the thermodynamic grand canonical potential w/T* = —{\nZ)/{VT^). Also, pressure, which 
is given by uj itself, is evaluated by integrating a derivative of w in current studies of the equation of state. Therefore, 
the calculations of the derivative of lj is basic for the study of QCD thermodynamics by lattice simulations, and the 
Taylor expansion method calculating higher order derivatives in fiq is the most natural extension from the study at 
fiq = to finite fiq. 



A. Taylor expansion of the grand canonical potential 



We study pressure p and quark number densities and defined by derivatives of the partition function 
Z{T, 

where fiu and fid are the chemical potentials for the u and d quarks. Let us define the quark chemical potential 
fjLq = {jjLu + lid)/'2- and the isospin chemical potential fii = — fid)/'^- Taylor expansion coefficients of physical 
quantities are given by derivatives of them in terms of fiu and Hd, or equivalently fiq and fij. We evaluate these 
coefficients at /i„ = /i^ = and study the physical quantities as functions of T and fiq in the isosymmetric case 
fJ-u = fJ-d = fiq (i-c. fii = 0). 



We define the susceptibility of quark number by 

Xq _ ( d d \nu + nd 



T2 \d{fi^/T) d(fid/T)J T3 ' 
and the susceptibility of isospin number by 

XI f 9 d \ Uu - Ud 



T2 \d{fiu/T) difidiT)) T3 



(9) 



(10) 



These susceptibilities correspond to the fluctuations of baryon number and isospin number in the medium, respectively 
[5^. They are expected to behave quite differently near the critical point in the (T, fiq) plane. 
We define the Taylor expansion coefficients of the pressure p{T,fiq) for the case fiu = fid = fiq as 



-2^c„(T)(-^l , c„(T) - 



(11) 



Here, co(T) is the pressure at fia = and has been computed by the CP-PACS Collaboration with the same action 
on 16^ X 4 and 16'^ x 6 lattices 0, j8|. Its value in the quenched limit is given in 51]. 

We also expand the quark number and isospin susceptibilities for the case fi^ = fid = fJ-q' 

XqiT,flq) _^ , (i^qV , Xl{T,flq)_^l 



= 2C2 



rp2 

where 

, _ 1 Ni 9" In Z(T, fiq + flj,flq- fij) 



I2c,{^)\-, ^i^ = 24 + 12cl(^)- + ..., (12) 



n\ Nf d{fii/TYd{fiq/Ty 



(13) 



1. Free quark-gluon gas at high temperature 



We expect QCD in the high temperature limit is described as free gas of quark and gluon. The pressure of the free 
jas in the continuum theory is given by 



P_ 

rp4 



87r2 
45 



E 

f=u,d 



7tt^ 1 /fifY 1 nif 

60 2 V T / 47r2 V T 



(14) 
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Note that the /Xg-dependence appears only through terms of fi'^ and /i^. The quark number density is a cubic function 
of Hq too. The quark number and isospin susceptibihties are the same for the free quark-gluon gas and are given by 
a quadratic function 



(15) 



Therefore, the Taylor expansion will converge well in the high temperature region. 



2. Hadron resonace gas at low temperature 



On the other hand, QCD at low temperature may be modeled by free gas of hadron resonances [53 ■ The partition 
function of the hadron resonance gas consists of mesonic and baryonic contributions, 



In Z(T, V, iJiq) = ^ In Z^^ (T, V, fig) + J] In Z^^ (T, V, 

i£ baryons 



2£ mesons 



where 



with energies ef — k'^ + mf and fugacities 



exp((3B,/i,)/T) . 



(16) 



(17) 



(18) 



Here Bi is the baryon number: Bi = 1,-1 and for baryons, anti-baryons and mesons, respectively. The upper sign 
in Eq. ([T7)) is for bosons, while the lower sign for fermions. Note that Z^, is actually independent of fiq. Expanding 
the logarithms in powers of fugacity, the integration over momenta, fc, can be carried out: 



27r2 



1=1 



(-1) 



i+i 



IrrLi 



(19) 



where K2 is a modified Bessel function. For rrii ^ T, the Bessel function can be approximated by K2{x) ^ 
^/tt/2x e~^(l + 15/80; + 0{x~'^)). Terms with £ > 2 in the series given in Eq. (flQl) thus are exponentially sup- 
pressed. 

Let us study the ^q-dependence of the partition function. The mesonic sector has no //^-dependence because Bi — 
for mesons. On the other hand, the baryonic sector can be approximated by the leading term in the expansion of Zi, 
since all baryons are heavier than a typical temperature scale. We obtain 



piT,fiq) p(r,o) _ 1 



ji4 



rp4 



[lnZ(r,M,)-lnZ(r,0)]~F(T) 



cosh ( ^ 



- 1 



with 



baryons 



T 



(20) 



(21) 



Note that each term in the sum for F now counts both baryons and anti-baryons. The quark number susceptibility 
is then given by 



Xq_ 
rp2 



9F{T) cosh 



T 



From Eq. (|20|) . the ratios of the expansion coefficients of p/T* in Hq/T are derived, 

C2n+2 9 



C2n 



(2n + 2){2n+l)' 



(22) 



(23) 



The ratio decreases as the order becomes higher. This means that the contribution from the higher order terms of 
/ig/r is small in the region of /ig/T^ 0(1). 
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FIG. 6: Time history of Di x (NfN^Nt)'^ (top) and V2 x {NfN^Nt)~^ (bottom) obtained by diflFerent noise sets at T/Tpc = 
0.925, mps/mv = 0.8. 



3. Numerical study near the transition temperature 



The behavior near the transition temperature is non-trivial. We expect a critical point at finite Hq. The Taylor 
expansion must break down at that point. We perform numerical simulations to study the expansion coefficients near 
the transition point. Using fi = Hqa, the explicit forms of the Taylor expansion coefficients are 



C2 



Nt 
2m 



A2, C4 



1 



2^ 



Nt 
2m 



B2, ci 



1 



■(S4- 62-42), 



(24) 



with 



I.e. 



A2 
B2 



{V2) + (Vl) , Aa - (2?4> + 4 {V3V1) + 3 (Vl) + 6 {V2VI) + (Vf) 
(2?2> , ^4 = (Vi) + 2 + (pi) + {V2V\) , 



a" IndetM 



/- 



Vi = Nf 



V2 - Nf 



12tr A/-i_M~i — M-i — 
\ nil. nil ni 



■T^- IT^ + 2tr M-i — M-i— M-i — 
Ofi a^x'^ J \ a/i Oyu ^ 



The derivative of the fermion matrix M at /i = is 



dfj. / \ dfi dfi J 



K ({1 - j,)U,{x) 5,+4,,-(l 
-^((l-74)(74(a;) 5^+4,, + (1 



74){/4^(a;-4) 
.+4,, + (l + 74)t/4^(x-4) 




for n : odd. 
for n : even. 



(25) 



(26) 



(27) 



(28) 



B. Random noise method 



We apply a random noise method to evaluate the traces in Eq. (|27l) . As we will see later, this method is effective 
when off-diagonal elements of the matrix are small. Therefore, the method works well for traces over spatial indices: 
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FIG. 7: Time history of the imaginary part (top) and real part (bottom) of Vi x {N f Nf Nt) ^ obtained by different noise sets 
at T/Tpc = 0.925, mps/mv = 0.8. 



Because the inverse of the quark matrix M~^{x, y) decreases as a function of |a; — yj, the off-diagonal elements in the 
spatial coordinate will be smaller than the diagonal ones. The random noise method will work well to suppress these 
small contaminations of off-diagonal elements. On the other hand, the off-diagonal elements in the color and spinor 
indices at the same spatial point are not suppressed by jx — y|, and will have the same magnitude as the diagonal 
elements. Because a staggered-type quark does not have the spinor index at a spatial point, the number of off-diagonal 
elements is only 6 in the 3x3 matrix, the contamination of off-diagonal elements may be not so serious. However, 
for Wilson-type quarks, because the number of color-spinor index is 3 x 4, the number of the off-diagonal elements 
in the quark matrix is 11 times larger than the diagonal one, so that the color-spinor index should be treated more 
carefully with Wilson-type quarks. In this study, we apply the random noise method for the spatial coordinates only, 
repeating the calculation for each of the color and spinor indices. 
We generate noise vectors (77^.0)^ ^ = ?y(«, x) (Sq.^, which satisfy 

^ 77(z,x)77* «<5,,, (29) 

-^^noisc 1 

for large iVnoisc- We adopt U(l) random numbers as 77, which are complex random numbers with jr/j = 1 and are 
generated from uniform random numbers 6 £ [0, 27r) with 77 — e*^. For each color-spinor index {a — 1, • ■ ■ , 12), we 

generate A^„oiso noise vectors (i = 1 A^noiso)- Then \imN„^.^,^oo{l/Nnoisc) Yliti"" Sali {'ni,a)^^p {vla)y ^ = ^x,y5i3,-i, 
hence 

gS'-iP^*- ("=1.V-). (30) 

where Xi^a = M~^(i9"^M/i9/i"^) • • • M~-^rii,a- To obtain X, we solve equations MXn — Yn recursively with Yi — rj, 
Y2 = {d''M/dn'')M-^ri = (9"M/9/x")Xi, etc. 

Because N^^.^^J2iVihx)T]*ihy) is O{y^l/Nnoisc) for x ^ y, errors due to finite iV„oisc decrease as ©(v^l/iVnoisc )• 
However, these errors are produced from all off-diagonal elements of the matrix in Eq. (I30p . hence these are proportional 
to the magnitude and number of the off-diagonal elements. Therefore, when the off-diagonal elements are not smaller 
than the diagonal elements, a number of noise vectors are needed to remove the error. This is the reason why we do 
not use the random noise method for the color-spinor index. 

For a product of traces, the random noise vectors for each trace must be independent. We compute such product 
by subtracting the contribution of the same noise vector from the naive product of two noise averages for each trace. 
This effectively increases the number of noises to iVnoiso (-/Vnoiso — 1) for the products and thus suppresses their errors 
due to the noise method. 

We then average over configurations to evaluate the expectation values in Eq. (^5)) . In addition to the errors due 
to the noise method, the statistical fluctuation of configurations contributes to the final error. To check the relative 
amount of the errors from the noise method, we calculate the operators I?„ (n = 1-4) using two independent sets of 
noise vectors with iVjjoiso = 10 on the same configurations. Figure [5] shows the time history of the imaginary part of Vi 
and the real part of 2?2 computed using these two sets of noise vectors. The operator I?„ is real for even n and purely 



12 



SB (^=41 



SB limit 



0.5 1 1.5 



J , I 1 I 1 I , u 

2 2.5 3 3.5 4 

T/T 

pc 



1 1 1 1 1 1 

* X IT^ H!pg/my=0.80 


SB (W =4)- 






- • • 

■ 

• 


• • 

- 


■ 
• 


SB limit 


■ 












1.5 2 

T/T 

pc 



FIG. 8: Quark number (black) and isospin (red) susceptibilities at /j, = for mps/my = 0.65 (left) and 0.80 (right). 
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FIG. 9: The second derivatives of quark number (black) and isospin (red) susceptibilities at /.tg = for mps/mv = 0.65 (left) 
and 0.80 (right). The dashed line is a prediction from the hadron resonance gas model: d^Xq/dl^ ~ Sxg/T^ 



imaginary for odd n [22l |. Therefore, the average of 2?i is zero because the expectation value is always real at /i^ = 0. 
We find that two results of 2?2 obtained by different noise sets are consistent with each other on each configuration, 
while two results of I?i are sensibly different. This means that, in the evaluation of I?i with A'noisc — 10, the error 
from the noise method is larger than the error from the statistical fluctuation of configurations. We can reduce the 
error from the noise method by increasing iVnoise- We plot the time history with A^noisc = 200 in Fig. [T] Two results 
of Im[2?i] using different noise sets are almost consistent, i.e., the error in I?i is now dominated by the statistical 
fluctuation of configurations with this A^noise- 

The required number of noise vectors depends on each operator. Here, we note that, in the evaluation of C4 and 
C4 through Eq. (|25p . the errors due to the error of I?i is dominant. In order to efficiently reduce the total errors of 
C4 and C4, we adopt large iVnoiso only for 2?i, keeping TVnoisc for other operators small. The values of A^noiso we adopt 
are summarized in Table HVl We choose A^noiso = 10 for the calculations of the operators in Eq. ([25]) except for the 
operators tr[(9"M/(}/i")Af ~^], where n = 1 — 4, for which we adopt A'noisc — 100-400 (the first number in the column 
of TVnoise in Table |IV|. 

Finally, we take advantage of the knowledge that the odd derivatives are purely imaginary and the even derivatives 
are real. In the lower panel of Fig. [71 we plot Re[I?i] which should vanishes when iVnoiso is large enough. We find 
that, unlike the case of Im[I?i] shown in the upper panel of the same figure, Re[2?i] data from two sets of random 
noises show no correlations in the time history even with small iV„oisc- Therefore, to further reduce errors from the 
random noise method, we can put the real and imaginary parts of the odd and even derivatives to zero, respectively. 



C. Quark number density, quark number susceptibility and isospin susceptibility 



We perform a series of simulations along LCP's for two quark masses corresponding to mps/rny — 0.65 and 0.80 
to calculate the expansion coefficients 02,04, C2 and c{ defined in Eq. ((24|) . The results are summarized in Table HVl 
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TABLE IV: Results of the Taylor expansion coefficients for mps/mv = 0.65 and 0.80. The first number in the column of A'l 
is A'noisc for the calculations of tr[(9"M/9/i")A'f^^], and the second number for other traces. See text for details. 





C2 X 2 


C4 X 4! 


ca X 2 


c\ X 4! 


iioisG 






mps/mv 


= 0.65 








0.82(3) 


0.352(59) 


6.3(108) 


1.189(6) 


1.41(49) 


400, 


10 


0.86(3) 


0.420(71) 


2.6(154) 


1.392(6) 


1.81(46) 


400, 


10 


0.94(3) 


0.963(64) 


10.5(103) 


1.857(10) 


2.88(64) 


400, 


10 


1.00(4) 


2.134(53) 


24.4(107) 


2.780(21) 


7.83(111) 


200, 


10 


1.07(4) 


4.140(27) 


8.7(21) 


4.396(16) 


5.58(34) 


200, 


10 


1.18(4) 


4.732(21) 


7.8(11) 


4.910(8) 


4.82(19) 


200, 


10 


1.32(5) 


4.938(20) 


7.1(14) 


5.052(6) 


4.65(13) 


100, 


10 


1.48(5) 


5.042(17) 


5.6(12) 


5.143(6) 


4.72(14) 


100, 


10 


1.67(6) 


5.133(15) 


4.0(11) 


5.229(5) 


4.67(13) 


100, 


10 


2.09(7) 


5.314(11) 


5.0(6) 


5.368(4) 


4.65(8) 


100, 


10 


2.59(9) 


5.447(13) 


4.8(6) 


r' A oc\ / a\ 

5.482(4) 


4.72(5) 


100, 


10 


■i.ZZ{lZ) 


5.517(12) 


6.4(7) 


5.562(4) 


5.05(8) 


100, 


10 


4.02(15) 


5.593(12) 


5.8(6) 


5.618(4) 


5.03(7) 


100, 


10 






mps/mv 


= 0.80 








0.76(4) 


0.066(34) 


3.8(51) 


0.549(4) 


0.37(19) 


400, 


10 


0.80(4) 


0.134(33) 


1.9(39) 


0.637(5) 


0.35(23) 


400, 


10 


0.84(4) 


0.251(35) 


0.0(37) 


0.776(6) 


0.80(27) 


400, 


10 


0.93(5) 


0.713(40) 


2.0(48) 


1.313(9) 


1.94(34) 


400, 


10 


0.99(5) 


2.071(34) 


17.4(47) 


2.498(17) 


5.13(53) 


400, 


10 


1.08(5) 


3.877(19) 


8.0(10) 
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4.92(18) 
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10 


1.20(6) 
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7.8(9) 


4.508(7) 


4.63(14) 
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10 
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10 
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5.8(3) 


5.234(5) 
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200, 


10 
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5.9(3) 
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10 


3.01(15) 
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10 



The results for Xql^"^ and xilT""^ at /ig = are plotted in Fig. |8l The circle and square symbols are for Xql^''^ ^^^^d 
XI /T"^, respectively. The short lines in the right end denote the values in the free quark-gluon gas (Stefan-Boltzmann) 
limit, both for Nt — A and in the continuum (cf. Appendix \K\ . 

At iiq = 0, Xq/T"^ = 2c2 and X//T^ — 2cl. Because Vi is a pure imaginary number, is negative in Eq. (1^51) and 
thus x//T^ will be larger than x^/T^, while the difference should vanishes in the high temperature limit according to 
Eq. (HH) for the free quark gluon gas. In the low temperature phase, Xq/T"^ and X//T^ correspond to the fluctuations 
of baryon and isospin numbers, respectively. Since the fluctuation of isospin number is mainly caused by pions, the 
fluctuation should be larger than that of the baryon number. Moreover, because the pion mass is more sensitive to 
the quark mass than baryon masses, Xi/T"^ "^iU show more sensitivity to the quark mass than x^/T^. 

As seen from Fig. [U both Xq/T"^ and xi/T"^ increase sharply at Tpc, in accordance with an expectation that the 
fluctuations in the quark-gluon plasma phase are much larger than those in the hadronic phase. We find that xi /T"^ 
is larger than Xq/T"^ at low temperatures and the difference vanishes in the high temperature region. Also, the isospin 
susceptibility increases as mps/nzy decreases at low temperatures, while Xq/T'^ does not ch ang e very much. These 
results agree quahtatively with previous results obtained with staggered- type quarks [TtI. [TsI. [20[[21|. [ssj . 

The quark number and isospin susceptibilities are expected to show quite different behaviors near the critical point 
at finite density. When the quark mass is nonzero, iso-triplet mesons are massive and thus are irrelevant to the critical 
behavior. Therefore, the iso-triplet susceptibility xi '^il^ not show singularity. On the other hand, if there is a critical 
point in the (T, /i^) plane, scalar sectors, il^ip and ipjo'4'i niay become massless at the critical point. We then expect 
divergence in the fluctuations of the chiral condensate and quark number towards the critical point. 

Figure [9] shows our results for d'^{xq/T^)/d{fiq/Ty\^^^^ = 24c4 (circles) and d'^ixi/T^)/difiq/T)^\^^^^ = 24ci 

(squares). We also plot dxq/T"^ as a dashed line in this figure to compare with the prediction from the hadron 
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FIG. 10: T-dependence of the pg-dependent contribution to the pressure, Ap/T* = p{^,q)/T'^ ~ p{0) /T'^ , at mps/mv ~ 0.65 
(left) and 0.80 (right). To is Tpc at ^i, = 0. 







FIG. 11: Quark number density at finite fiq for mps/my ~ 0.65 (left) and 0.80 (right). To is Tpc at /iq — 0. 



resonance gas model in Eq. (22), i.e. d^Xq/9fJ-q ~ ^Xg/T^- These results are consistent within the error at T < Tpc . 

Although the statistical errors are not quite small yet, the two susceptibilities show quite different behaviors near 
Tpc- d'^{xq/T'^)/d{fj,q/T)'^ near Tpc is more than three times larger than that at high temperatures, suggesting the 
large enhancement in the quark number fluctuations as the density is increased. Moreover, the peak height is larger 
for smaller mpg/my. On the other hand, no such sharp peak appears for d'^{xi/T^)/d{iiq/T)^, in accordance with 
the expectation that xi is analytic at the critical point. These observations suggest the existence of the critical point. 
Similar results were obtained by p4- improved staggered fermions [13, 0, [13, El . 

Finally, we evaluate the equation of state at finite /ig combining the results of derivatives. Figure [TU] shows the 
yitq-dependent contribution of the pressure, Ap/T'^ = p{fiq)/T'^—p{0)/T'^ = C2{fJ.q/Ty + C4{fj,q/T)'^, at mps/m-v = 0.65 
(left) and 0.80 (right). The truncation error is 0(/i^). Tq is Tpc at fiq = 0. The finite density correction for p/T"^ 
becomes the same size as p/T"* at /ig = around piq/T ^ 0(1), and the correction ApjT^ increases rapidly around 
Tpc in comparison with the behavior oi pjT^ at /ig = 0. This suggests that the pressure changes more sharply as 
[Xq is increased. The quark number density, nqjT'^ = (n^ + nd)/T'^ — 2c2(/ig/T) + 4c4{^q/T)^ + 0{p.q), is shown in 
Fig. 111! The quark number susceptibility and isospin susceptibility are shown in Fig. 1121 and Fig. I13[ respectively. As 
discussed above, we find large quark number fiuctuations near Tpc when /ig is increased. On the other hand, such an 
enhancement around Tpc is not visible in the isospin fiuctuations. These results are consistent with the observations 
with staggered- type quarks and suggest a critical point at finite fj,q. 



IV. EQUATION OF STATE FROM THE GAUSSIAN APPROXIMATION 



In the previous section, we have studied the equation of state at finite density by computing the Taylor expansion 
coefficients c„ up to the fourth order, based on the calculation of I?„ = iV/[9" Indet for n < 4. We found, 

however, that the statistical errors in Uq/T^ and Xq/T"^ ^-re not small. Furthermore, the statistical errors will be 
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FIG. 12: Quark number susceptibility at finite /i, for mps/mv = 0.65 (left) and 0.80 (right). 
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FIG. 13: Isospin susceptibility at finite /i, for mps/my = 0.65 (left) and 0.80 (right). 



larger when we include higher order terms, Cg, Cg etc. 

In this connection, we recall that, in a previous study with staggered- type quarks j23j . a hybrid method ol the 
reweighting technique and Taylor expansion 22], combined with a Gaussian approximation for the complex phase 
distribution of quark determinant, was efficient to suppress statistical fluctuations at finite densities. We call the 
method simply the Gaussian approximation. In this section, we apply the Gaussian approximation to the calculation 
of EOS with improved Wilson quarks. 

In the evaluation of higher order Taylor coefficients c„ with n > 4, the calculation of 2?„ at large n is quite 
demanding. However, the free quark-gluon gas leads to 2?n = for n > 4 in the continuum limit. Therefore, we may 
approximately evaluate higher order coefficients by keeping 2?„ for ti < 4 only. The approximation should work at 
least at high temperatures. Therefore, we consider the following approximate grand canonical potential. 



rp4 



1 



1 



In 



VU (detM(^)) 



„-s. 



lnZ(r,0) 



In ( exp 




lnZ(r,0) 



(M=0) 



1 ^^// detM(M) \^' 
VT^ ^ \ ^etM(O) J 



(m=o) 



(31) 



where ^ = fiqU = fiq/(TNt) and iV,nax = 4. Here, 



/(m=o) 



is the average over configurations at /i = 0. This 



approximate grand canonical potential is equal to the exact potential up to 0(//^°""=), and most of higher order 
contributions are contained except for terms including I?„ for n > A^max- In this context, the method would be better 
than the truncated Taylor expansion method discussed in the previous section. 
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FIG. 14: The histogram of Q for simulations at {m^s/mY ,T /Tj,^) = (0.65,0.94) (top left), (0.65, 1.32) (top right), (0.80,0.93) 
(bottom left) and (0.80, 1.35) (bottom right). 
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FIG. 15: The distribution in the {F,e) plane for pq/T = 0.5 (left) and 1.0 (right) at (mps/my, T/Tpc) = (0.80,0.93). 



A. Gaussian approximation for the 9 distribution 

We calculate the grand canonical potential (|3T|) following the method of Ref. [IS]- We first rewrite the grand 
canonical partition function as follows. 

2(T.,., = Z(T,0)/(|iii<|)"'\ .Z(T.O,(.™.*..) (32, 

where F{^) and 0(/Lt) are the real and imaginary parts of A^f ln(detM(//)/detM(0)), respectively, and they can be 
calculated by the Taylor expansion in /i. Since odd (even) derivatives of ln(det M(/i)/detM(0)) are purely imaginary 
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FIG. 16: {e\F - {F)))/{{F - {F)f) for fig/T = 0.5 and 1.0. 



(real), we have 



F(^) = iVfRc 



/ detM(M) \ 
^detM(O) / 



1 



ra2"(indetM) 



{ (2n) 



■Re 



9^' 



2ri 



,2n 



oo 



(m=o) 



^ (2n) 

n— 1 ^ ^ 



■ReI?2nM 



2n 



(33) 



In this paper, we study terms up to /i'*. For the complex phase 9, we have 



0{fi) = 7VfIm[lndetM(/i)] 



1 



[2n+l) 



■Im 



a2"+i(lndetM(^)) 



2ri+l 



,2n+l _ 



(m=o) 



E 



„t^o(2n+l)! 



ImI?2n+lAi 



2n+l 



(34) 



We note that lndetM(/i) is not uniquely defined for complex det M{fi). On the other hand, the /i derivatives of 
lndetM(/i) are unique. We regard the Taylor expansion in (p4)) as our definition of 9. Note that the 9 thus defined 
is iVOT restricted to be in the range — tt to tt, and the maximum value of \9\ is infinite in the large volume limit. The 
principal value of N{ In det M[^) is recovered by identifying 9 + 2mT with 9 in the range — tt to tt. 

Histograms of 6* are shown in Fig.[Tlfor /i^/T = 0.5 and 1.0 at (mps/my, T/Tp^) = (0.65, 0.94) (top left), (0.65, 1.32) 
(top right), (0.80, 0.93) (bottom left) and (0.80, 1.35) (bottom right). We find that the fluctuations in 9 become larger 
as /iq increases. Note that the width of the distribution is larger than 27r at T < Tpc- A large fluctuation in 9 makes the 
calculation of In Z (T, /i^ ) difficult due to a rapid change of the factor e*^ . This is the origin of the sign problem. On the 
other hand, these figures suggest that the distribution of 9 defined in this way is almost of Gaussian. In Sec lIV Bl we 
discuss that the Gaussian approximation corresponds to the leading order approximation of the cumulant expansion 
and confirm the validity of the Gaussian approximation. This is a key observation to avoid the sign problem: In a 
previous study with staggered quarks, using the fact that the ^-distribution is well described by a Gaussian form, 
the ^-averaging has been carried out. The resulting errors for observables turn out to be smaller than those with the 
naive averaging, and thus the method may enable us to perform a reliable evaluation at a wider range of [23| . 

To implement this assumption, we define the distribution function w{F, 9) as 



w{F, 9) EE VU S{F - F{n)) S{9 - 9{p)) [det M{0)] 



Z{T,Q){5{F ~ F{y)) 5(9 - 9(p)))(^,^o) 



(35) 



where 9{ii) and F{^) are defined in Eq. (|34l) and Eq. (|33p . Note that w{F^ 9) depend implicitly on /i. Figure [15] shows 
a typical distribution of {F,9) at {mps/^v,T/Tpc) = (0.80,0.93). The Gaussian ^-distribution means that 



^u{F,9)^^l^^woiF)e-^^-^^^'' 
With this form, it is easy to carry out the 6'-integration as follows. 



(36) 



dF I d9 w{F, 9) e^e^" 



dF I d9 e^wa{F) 



a2{F) 



= / df^e^i«o(F)e-i/(4''^(^» = Z(T,0)(e^('')e^i/(4'^^(^(^»\ 



(37) 
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In the last line we use the fact that 



wo{F) = / VUS{F~Fiti)) [dctM(0)]^^'-e-^« =Z(r,0)((5(F-F(/i)))(^^o) 



(38) 



holds within this assumption. Note that the problematic factor e*^ in Eq. ([5^ is now replaced by a positive definite 
factor e~^/('*"=^). Thus the statistical error of Eq. (j37p is always smaller than its expectation value, i.e. there is no 
sign problem. 

Of course, one may replace the Gaussian distribution function w{F,9) with a periodic distribution function given 

by 



1 ^ 

lim — V w(F,e + 2TTn). 

n=-N 



(39) 



However, the integral of e*^ does not change simply because / e'^^w{F,9 + 2TTn)d9 gives the same answer as 
J e^^w{F,d)d9. Hence, the absence of the periodicity of 27r in w{F,d) is not a problem for the integral of e*^. 

The validity of this method can be discussed more precisely based on the Taylor expansion of the partition function 
at least in the low density region. In Appendix |B] we compare the derivatives of In Z in the Gaussian approximation 
with the exact calculations up to 0(/i^). We find that the Gaussian approximation does not affect up to 0(/i^). At 

the fourth order in /iq, {"Df) of Eq. ([27| is replaced by 3 {'D'l)^ in the Gaussian case. In Ref. [HJ], the effects caused by 
deviations from the Gaussian distribution in w{F,6) are estimated assuming w{F,9) ^ exp[— a20^ — a^O'^]. It turned 
out that the additional term 04 does not affect the terms up to /i^ as far as ai/a2 < 0{1). 
Now the problem is reduced to a determination of the coefficient 02 {F) : 



2a2(F) 



(g2(/i)J(F-F(^)))^^^g^ ^ /Pt/ §{F-F{^l)){deiM{Q))^n 
{6iF-Fi^))) 



(m=o) 



JVU (S(F-F(^))(dctAf(0))^f£ 



(40) 



The distribution shown in Fig. [15] suggests that the i^-dependence in (0^) p is mild. Unfortunately, the limitation 
of the statistics makes a precise evaluation of {0'^)f for each thin slices of F difficult. However, when we restrict 
ourselves to calculate the equation of state up to 0(/x^), we only need to evaluate the first derivative of (0^) f in terms 
of F: Because I?i and I?2 represent the leading /ig-dependence of 9 and F, respectively, consulting Eq. (|25|) . we note 
that the i^-dependence of {9'^)f affects only in the {T>2'D\) term for the 0(/iq) coefficients C4 and c\. (See Appendix 
|B]too.) This quantity, i.e. the 0(/i^) contribution of {F9'^), corresponds to the first derivative of {9'^)f because 



'■{^^)(F{^,)-{F)))^,^o) = 



')iF){F-{F)) + 



dF 



{F-{F)r + 



wojF) 
Z{T,0) 



dF 



d{9^ 



dF 



{{F^{F)r)^,=o), 



when the i^-dependence in 
respect to F as 



) F is mild. Using this relation, we then estimate the first derivative of 



d{9^ 



dF 



(F) 



{0'iF-{F))) 

{{F-{F)y) ■ 



(41) 



with 



(42) 



which is shown in Fig. [TOl We find that [d{9'^) f / dF](^F) is actually smaller than statistical errors, so that 
(^^)(F) is a good approximation. This point is also suggested in chiral perturbation theory |55i] . To include the small 
F dependence of 02 (-F), we assume a simple ansatz function; 



1 



2a2(F) 



)i=^ = /(^) = cxp[a;i+a;2F], 



(43) 



where we take into account the fact that 9'^ is positive for all F. The two parameters are sufficient for the exact 
calculation up to 0(/i^). We thus determine fit parameters xi and X2, by minimizing = '^j\9f — /(i^i)]^, where 
the summation is taken over configurations. 

Finally, we integrate over F. The factor in Eq. ([57)1 is a potential danger in the integration because it can easily 
shift the central contribution for the average to a statistically poor region of F . This will be the case when /ig is 
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not small ((F) is not small). At small ^, this problem can be removed in part by a reweighting in the /3-direction of 
the coupling parameter space such that the fluctuation in e^'-'^-' is compensated by that in the gauge action. This is 
possible since F is strongly correlated with P = — S'g/(6iVsitc/3), where the gauge action Sg is defined in Eq. ([3]), and 
A'sitc = X Nt- By reweighting, the expectation value of an operator O at /? is evaluated from a simulation at /3o as 



/p6Ar3itc(/3-^io)P\ 
\ ho 



- /^fiAf_,._r«_fl„^P\ ~- (44) 



To calculate (e^('")e^i/(''°"(^))), we adjust j3 such that the value of e^e-i/('*°")e'^^-"=('^-'3«)^ is stabilized during the 
Monte Carlo steps. In practice, since e^'^^'^'e"^/'*"^'-''^'' — 1 bX — we start with /3 = /3o at /i^ = and find /3 for 
finite /i, at which the fluctuation of e^e-i/('*°2)e6^-'=('3-'3«)^ = X, 

:^-W(.=o))')^^^^y WU), (45) 

is minimized. Since F becomes larger for larger P , (3 < The resulting shift in /3 is translated to the temperature 
scale using a cubic spline interpolation of the temperature data. Because we do not shift the hopping parameter, a 
shift in /? leads to a slight deviation from the original line of constant physics (LCP). In our study, however, the shifts 
in (3 turn out to be smaller than 0.03. Since these shifts are negligible in Fig. [1] we disregard the resulting small 
deviation from the LCP, and simply translate the shifts in /3 to shifts in T for the final plots. 
To conclude we summarize the final formulae: 



— , (fc* )i^ = exp(xi + a;2i^). (46) 



Z(r,0) (^^&N,,UP-MP) 



B. Gaussian approximation as the lowest order approximation of cumulant expansion 

The only difference between the Gaussian approximation (|46p and its exact formula is the replacement of (cxp(i0)) p 
by exp[— (^^)i?/2]. The meaning of the replacement can be understood in the context of the cumulant expansion. 



(exp iti) p = exp 



(47) 



where (6'")c is the n^^ order cumulant, e.g. 

{0'). = {e')p , {e')^ = {0'), - 3 {0% , (9^)^ =. (e^)^ - 15 (e^)^ {e^)^ + 30 {0% . (48) 



Note that (6'")c = for odd n due to the symmetry under 9 — > —9. Because only the odd-order cumulants are the 
source of the complex phase in (exp(i0))i?, the value of {exTp{i9)) p is guaranteed to be real and positive from this 
symmetry if the cumulant expansion converges. There is thus no source of the sign problem once we eliminate the 
odd terms. 

When the distribution of 9 is of Gaussian, the 0(9'^) terms vanish for n > 2 in Eq. (|47)) . Hence, the Gaussian 
approximation is equivalent to the approximation that the higher order cumulants are neglected except for the first 
nonzero term. If one wants to improve the Gaussian approximation, it is achieved by adding higher order terms. 

Moreover, the cumulant expansion can be regarded as a power expansion in terms of because 9 ^ 0{^q). 
Therefore, if we take into account the cumulants up to the n*^ order, the truncation error does not affect the Taylor 
expansion up to 0(/i^'). The Gaussian approximation corresponds to the leading non-trivial order approximation of 
the Taylor expansion in /ig. 

On the other hand, a careful discussion about the infinite volume {V) limit is required. Because the operator 9 is 
roughly proportional to V , the n*^ order cumulant (0")c may increase as 0(T^") naively. In such a case, the cumulant 
expansion does not converge at large V . However, the following argument suggests that the convergence property of 
the cumulant expansion is independent of the volume when the correlation length of the system is finite. Note that, 
since no critical point is expected to exist in two-flavor QCD at > and /ig = 0, the correlation length between 
quarks is finite. 

The expansion coefficients of 9 in Eq. ([51)1 are given by combinations of traces of products of Af^^, d'^M/d{iiq/T)'^ 
and so on. For example, 2?i is given by the trace of Nf[AI~^{dM/d{iJ,q/T))] and the diagonal element of this matrix 
is the local quark number density operator (~ i}}^Qip{x)) at pLq = 0. If the correlation length of the local number 
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FIG. 18: The relative magnitude of the fourth order cumulant contribution to the leading order contribution as a function of 
the temperature for mps/mv = 0.65 (left) and 0.80 (right). 



density operator is much shorter than the system size, we may decompose I?i into independent contributions from 
spatially separated regions. The same discussion can be applied to higher order coefficients I?„ too. 

In this case, one can write the phase as 9 = J^x where 0x is the contribution from a spatial region labeled by x 
and these contributions are independent. The average of exp(i6') is thus 

(e-)«n(^"^)--pfEES(^")c)- (49) 

X \ X n / 

This equation suggests that all cumulants {0")c ~ J2x (^")c increase in proportion to the volume as the volume 
increases. Therefore, while the width of the distribution, i.e. the phase fluctuation, increases in proportion to the 
volume, the ratios of the cumulants are independent of the volume. The higher order terms in the cumulant expansion 
are well under control in the large volume limit. 

Because 6 is 0{^q) and {0")c is 0(/i^), the Gaussian approximation is valid at small fig and the higher order 
cumulants will become visible at large fiq. The application range of the Gaussian approximation in terms of /ig 
must be checked for each analysis by calculating the ratio of cumulants. However, the volume-dependence of the 
ratios of cumulants suggests that the application range does not change once the system size becomes larger than 
the correlation length. This means that the qualification of the Gaussian distribution on a small lattice is enough to 
verify the Gaussian approximation. 

We study the validity of the Gaussian approximation by examining the relative magnitude of the fourth order 
cumulant contribution to the leading order contribution in Eq. (j47|) : 

The Gaussian approximation is valid if i? ^ 0(1) is satisfied. In this paper, we will check whether R is consistent 
with zero, which is a less stringent condition when the statistical error is large. 
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FIG. 19: /iq-dependent contribution to pressure as a function of T/Tq for each \iqlT with mps/m-v — 0.65 (left) and 0.1 
(right), where To is Tpc at Hq = 0. 



Here, we note a caveat in the evaluation of (^?^) from the histogram. Because we calculate 9 using the random noise 
method, the fluctuation of 9 contains a contribution due to the finite number of noise vectors (A^noisc)- This makes 
the width of the histogram wider than that of the true distribution. True width is given by in the limit of 

large A'^noisc- To reduce the errors in {9"^) due to finite iVnoise, we adopt the subtraction method discussed in Sec. HIT Bl 
for the calculation of products of traces. The expectation value of 9"^ is summarized in Fig. [iTl Filled symbols in 
Fig. [T7]are the results of the subtraction method. We have checked that the A'noisc-dependence in {9"^) is negligible 
with our choices of A^noisc- We find that (6'^) becomes larger than 0(7r^) from /ig/T ^ 0.5 in the low temperature 
phase while, in the high temperature phase, the complex phase fluctuations decrease as T increases, in accordance 
with our expectation that the quark determinant is real in the high temperature limit. On the other hand, the width 
of the histogram shown in Fig. 1 141 corresponds to (9^) obtained by the naive calculation without subtraction, which 
is plotted with open symbols in Fig. [171 The difference between the results by the subtraction and naive methods 
decreases as iVnoisc increases but is almost the same size for all temperatures, and the error due to finite A^noisc is larger 
than the expectation value of (0^) at high temperature. Therefore, the subtraction is indispensable for a calculation 
of the width of the 9 distribution. 

We summarize the results for R in Fig. [THl The circle and square symbols are the results for Hq/T = 0.5 and 1.0, 
respectively. Filled symbols are the results of the subtraction method, while open symbols are the results of naive 
calculations without the subtraction. Although errors become gradually larger as f^q/T increases and are as large as 
0(1) for Hq/T = 1.0, the central values of R are consistent with zero for all temperatures and f^q/T ^ . However, 
we need higher statistics to identify the actual magnitude of R and to check the validity range of the Gaussian 
approximation in terms of ^.q/T, which is left for future investigations. 



C. Results for the equation of state and quark number susceptibility 

In Fig. [191 we show the results for the /iq-dependent contribution to the pressure, Ap/T** — p{iiq)/T'^ — p(0)/T'*, 
obtained by the Gaussian approximation. Comparing with Fig. IIO) improvement towards larger /ig is clearly seen. 
We calculate the quark number density Uq and its susceptibility Xq by the following numerical differentiations: 

nq _ Nf d{\nZ) Xq _ N! d^jlnZ) 

r3 Nid{f,q/Ty N!d{fiq/T)^- ^ > 

Results of ln[Z(r, fiq)/Z{T, 0)] around representative points fiq/T — 0.2, 0.4, • • • , 1.2 are shown in Fig. [20l where (3 
is optimized at each fiq/T. The value of \n[Z{T, ^q)/ Z{T,0)] increases as T/Tq increases for each fJ,q/T, where Tq is 
Tpc at fiq/T = 0. In Fig. [20l results at the optimized values of T/Tq{P) for simulations listed in Table U are shown. 



^ Because the complex phase vanishes in the high temperature hmit, [8'^) becomes smaller as T increases. The small {9^) causes the large 
statistical error of R at large T for the subtraction method. Where (6^) is small, however, the correction due to the phase fluctuation 
itself is small, and thus a deviation from the Gaussian approximation does not affect the results. 
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FIG. 20: /Xq-dependence of ln[Z{l3, ^q)/Z{l3,0)] for each temperature. The values of \n[Z{P, iJ.q)/Z{l3,0)] increases as T/To 
increases for each i-iq/T. The left and right figures are the results at mps/mv = 0.65 and 0.80, respectively. 




FIG. 21: Quark number density for each fiq/T at mps/mv ~ 0.65 (left) and 0.80 (right). To is Tpc at fiq = 0. 



We then fit the data in the range flq/T — 0.05 < Hq/T < jiq/T + 0.05 by a quadratic function of ^uqjT , 



Nl 



In 



Z{T,0) 



y3 y 



2T2 



C(m, 



(52) 



with the fit parameters nq{jlq),Xq{i^q) ^^id C{jlq), for each values of fiq/T and T/Tq. 

The results of nq{fj.q) and Xgif^g) £^re plotted in Figs. [2T] and [22l As is similar to the case of p/T^, the statistical 
errors in these figures are much smaller than the results given in Sec. IIIII Moreover, although simulations at different 
temperature are independent, the temperature dependence in these figures is smooth and natural. The reduced 
statistical fluctuations over the results of Sect. IIII Cl are mainly due to the Gaussian method for the 6'-averaging and 
the /3-reweighting for the i^-averaging. 

At TOpg/mv = 0.65, we find a sharp peak in Xg/T^ near Tpc- The peak becomes higher as fiq increases. These 
observations are consistent with the findings in Sec. IIIII and suggests a critical point at finite ^q. On the other hand, 
the peak is much milder at mps/my = 0.80. This may be explained in part by the expectation that the critical 
point locates at larger fiq because the quark mass is larger than that for mps/my = 0.65. Further studies with 
increased statistics around Tq are needed for a more definite conclusion. A scaling analysis increasing the volume is 
also important. 



V. HEAVY-QUARK FREE ENERGY AND DEBYE SCREENING MASS AT FINITE TEMPERATURE 

AND DENSITY 

In this section, we investigate the heavy-quark free energies between static quark (Q) and antiquark (Q), and 
between Q and Q. These free energies arc important inputs in phenomenologial studies of color-singlet quarkoniums 
such as charmoniums and bottomoniums in QGP [56, 57J and of color non-singlet quark-quark states in QGP [5^. 
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' 

FIG. 22: Quark number susceptibility for each ^iq/T at mps/my = 0.65 (left) and 0.80 (right). 



Lattice simulations for QQ and QQ free energies in different color channels at = have been performed in Nf = 2 
QCD with the staggered fermion [sl, H^l and with the Wilson fermion [13, [Ml- In these works, Coulomb gauge fixing 
is employed to define the Polyakov-loop correlations in different color channels. Furthermore, the QQ free eiiergy at 
finite /ig has been studied with the staggered fermion by the reweighting method in the parameter plane [62] and 
by the Taylor expansion method [63j . Screening masses at finite have been also studied in dimensionally reduced 
effective field theory at high temperature [64|. 

Here, we extend our previous study with two fiavors of improved Wilson quarks at /ig = [2^ to finite /ig using 
the Taylor expansion method. Under Coulomb gauge fixing, we calculate the expansion coefficients of the heavy- 
quark free energies up to the 2nd order with respect to Hq/T for color-singlet QQ channel, color-octet QQ channel, 
color-sextet QQ channel and color-antitriplet QQ channel. The effective running coupling and Debye screening mass 
are also extracted by fitting the screened Coulomb form expanded as a power series of /ig/T, and compare with a 
prediction of the thermal perturbation theory. 



A. Taylor expansion of heavy-quark free energy 



The expectation value of an observable O for fiu = fJ-d = fJ-q is defined as 

(0>M, = I -DUO [det M{^^)f^ e-^^ (53) 

where fi = LL^a. For O which does not depend on fj,g explicitly, {O) can be expanded as a power series of /i = /igO 
as follows [63i|: The quark determinant is expanded as, 

[detAf(/i)]'^^ = [detM(0)]'^^(l + MiM + M2Ai'+ 0(^*3)), (54) 

with the expansion coefficients Mi = Vi, M2 = 5(^1+ ^2) , etc. using 2?„ defined by (|26p . Then, using the fact that 
the system is symmetric under /^g — > — /^g, (C)/j, can be expanded as 



l + (A/2>o/i2 ^ ' ^ ' 

= {O)o [1 + Oia*+(-(A/2)o + 02)m' + 0(a*')] , 

where (0)0 = (0)^^=0 and Oi is defined by 

O- (56) 
- {O)o ■ ^^^^ 

The heavy-quark free energies are defined by correlation functions between Polyakov loops, n{x) = Jl^^Ii U4{t,x). 
At a fixed gauge, the QQ correlation function can be decomposed into color singlet (1) and color octet (8) channels. 
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while the QQ correlation function into color antitriplet (3*) and color sextet (6) channels as follows |65l l66j: 

n^{r) - itrl]t(x)r!(y), (57) 

Q^{r) = itrl]t(x)trr!(y)-^trr!t(x)r!(y), (58) 

/2«(r) = ^trl](x)trl](y) + ^trr!(x)r!(y), (59) 

Q^'{r) - itrO(x)trO(y) - itrf](x)f](y), (60) 

where r = |x — y|. The free energy for color channel R {R = 1, 8, 6, 3*) is defined as 

e-^Hr,T,,q)/T^(^^R^^^^ (61) 

Above Tpc, we introduce normalized free energies {V^ ) by dividing the right-hand side of ((6T|) by 
(L)^^(i)*^ for QQ free energies and {L)'j^^ for QQ free energies, where L = tril. vanishes at r — oo. The Taylor 

expansion of with respect to /ig/T is given by 

V^{r, T, M,) = + (^) +v^{^y + 0{^^% (62) 

where 



= 0, ,64, 



for color singlet and octet QQ channels, and 



T \ I 



= -lM^h (66) 



T - 7v/^i +^4' ^^^^ 
^^^((M.)o4«)^-r2.«)+2|-|, (68) 

for color sextet and antitriplet QQ channels. Here fl^ — (J?^Af„)o/(J?^)o, and the tn is an n-th order coefficient of 
the Taylor expansion of the Polyakov loop: 

W., = 4 + ^i(^)+^2(^)' + 0(M^). (69) 

Note that the color singlet and octet channels do not have the odd orders in the Taylor expansion since the free 
energies for both channels are symmetric under ^.q — — /ig, i.e., the QQ free energies are invariant under the charge 
conjugation. 



B. Results for expansion coefficients of normalized free energies 



Heavy quark free energies are calculated in the high temperature phase on the lines of constant physics at tops /tov = 
0.65 and 0.80 (see Table Observables are measured every ten trajectories at each quark mass and temperature, 
and the statistical errors are estimated by a jackknife method with the bin size of 100 trajectories. 

The results for the expansion coefficients of the normalized free energies at m-ps/my = 0.65 are shown in Fig. [23] 
for the color singlet and octet QQ channels, and in Figs. [51] and [55] for the color sextet and antitriplet QQ channels. 
Those obtained at tops/tov = 0.80 are shown in Figs. 
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FIG. 23: Vq" (left) and V2 (right) for color-singlet and octet QQ channels above Tpc at mps/my = 0.65. 
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FIG. 24: «^ (left) and nf* (right) for color-sextet and antitriplet QQ channels above Tpc at mps/mv = 0.65. 



The ti(f 's shown in Figs.[23l[24l[26]and[27]are the normalized free energies at jiq — 0. The fact that, with increasing 
the distance r, and increase while Wg and Wg decrease represents the finding of our previous study "24"] that, at 
/ig = 0, the inter-quark interaction is "attractive" in the color singlet and antitriplet channels and is "repulsive" in 
the color octet and sextet channels. 

From these Figures, we note that, both around Tpc and at higher temperatures, the sign of is the same with 
that of Vq, whereas the sign of a is the opposite of that of v^: 

Vi ■ Vq > (only for QQ free energies), (70) 
v^ -vl^ < 0. (71) 

Because is absent for QQ free energies, this means that, in the leading-order of /i^, the inter-quark interaction 
between Q and Q becomes weak at finite /i^, while that between Q and Q becomes strong. In other words, QQ (QQ) 
free energies are screened (anti-screened) by the internal quarks induced at finite /iq. 



C. Screening properties at finite T and fiq 

At /iq = 0, the color channel dependence in the free energies was shown to be well absorbed in the kinematical 
Casimir factor at high temperatures 24] , as first noticed in quenched studies [gtI . [ssj . Therefore, we fit the normalized 
free energies by a screened Coulomb form. 



(72) 
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FIG. 25: V2 for color-sextet and antitriplet QQ channels above Tpc at mps/mv ~ 0.65. 
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FIG. 26: The same figures as Fig. [23] at mps/my = 0.80. 



where the Casimir factors = (X]a=i ' ^2)^ fo^' various color channels are given by 

C. = 4 C- = l C.i C»- = 4 (73) 

At small /iq, the effective running coupling acff(r, /ig) and the Debye screening mass mo{T, iiq) are expanded by 
powers of Hq/T: 



acS 



«0 + «i(^)+a2(^)' + O(/.^), (74) 
mo = rriDfi + mD,2 (^^^ +0(/i'^), (75) 

where we use the fact that the Debye screening mass does not have the odd powers in the Taylor expansion because 
it corresponds to the self-energy of the two-point correlation of the gauge field which is symmetric under fiq — > —fiq. 
Properties of ao(T) and mD,n{T) are discussed in [24l|. 

Expanding ([72]) with respect to ^J.q/T using ([74)1 and (|75)) . and comparing with the expansion (|62|) of the normalized 
free energies, we obtain the following relations: 

vo{r,T) = c^.^£(Zle-™i.,o(T)r.^ (76) 
r 

vi{r,T) ai{T) 



vo{r,T) ao{T) 
V2{r,T) ^ a^jT) 
vo{r,T) ao{T) 



(only for QQ free energies), (77) 
mD,2{T)r. (78) 
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FIG. 27: The same figures as Fig.[24]at mps/my = 0.80. 




FIG. 28: The same figures as Fig. [25] at mps/my = 0.80. 



Therefore, the expansion coefficients of aog and mo for each T can be calculated by fitting the normalized free 
energies for appropriate ranges of r. We chose the fit ranges to be 0.5 < rT < 1.0 for Eq. (1771) and 0.25 < rT < 1.0 
for Eq. ([75)) . In Appendix [Cl we study the fit range dependence of the fits, and find that the magnitude of systematic 
errors in the expansion coefficients due to the fit range are at most comparable with that of the statistical errors at 
T > 1.2Tp,. 

The results for the first order coefficients ai(T), which appear only for the color sextet and antitriplet QQ channels, 
are shown in Fig.l^Hfor TOpg/my — 0.65 (left) and 0.80 (right). The second order coefficients a2{T) and mn,2{T) are 
shown in Figs . [501 and [5l1 at mps/my ~ 0.65 (left) and 0.80 (right), respectively. Numerical values of these coefficients 
are summarized in Appendix [Cl 

From these Figures, we find that there is no significant channel dependence in these coefficients at high temperatures 
(T ^ 2Tpc), similar to the case of aQ{T) and mo^oiT) studied in [24!|. We note that mD,2{T) is positive at T l.bTpc 
which means that magnitude of the Debye mass becomes larger at finite densities in the leading-order of /ig. This is 
qualitatively consistent with results calculated with an improved staggered quark action for the color-singlet channel 
[63] ■ We also find that, although ai{T) remains finite even at T ~ 4Tpc, the magnitude of a2{T) is almost zero for all 
color channels at T <; 1.5Tpc. Therefore, to reduce statistical fluctuations in mo,2iT), we may assume a2{T) = in 
the fit ([75)) . The results are shown in Fig. [321 Smallness of the color-channel dependence became clearer. Numerical 
values for ai{T), a2{T), mD,2{T) and TO£i_2(T;a2 = 0) are summarized in Tables [VTlIXI together with x^/A^df for 
each fit. 
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FIG. 29: Qi(T) at mps/mv = 0.65 (left) and 0.80 (right). 
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FIG. 30: 02 (T) at mps/mv = 0.65 (left) and 0.80 (right). 
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FIG. 31: mD,2(T) at mps/my = 0.65 (left) and 0.80 (right). 
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FIG. 32: Results of moMT) assuming a2{T) = 0, at mps/mv = 0.65 (left) and 0.80 (right). 



D. Comparison with the thermal perturbation theory 



The 2-loop running coupling is given by, 

g,f(..) = /3oln(J)%|lnln(J)', (79) 

where k and A are the renormalization point and the QCD scale parameter, respectively. In the thermal perturbation 
theory the argument in the logarithms can be decomposed as k/A = {K/T){T/Tpc){Tpc/ A) where we adopt A = 

A^"^ ~ 261 MeV l0\ and Tp^ ~ 171 MeV [7]. We assume that the renormalization point k is in the range k = ttT 
to 3ttT. Therefore, 1721 can be viewed as a function of T/Tp^. In the leading order of the thermal perturbation theory, 
the Debye screening mass with 521 is given by 



1 /2 

m^r?{T,^^,) = .g2i(«) \ { I + ^) + ^M^j • (80) 



Thus, the leading-order expansion coefhcients are given by 



Taking the ratio of these coefficients we find for N 



f 



-.LO 
'■D,2 



3 



mLO, 87r2 ■ 



(82) 



In Ref. [24| . we found that, at /i^ — 0, the leading-order thermal perturbation theory predicts much smaller values 
for mD,o{T) than the lattice results. In the left panel of Fig. [331 '^e compare our results of to_d,2(T) for the color 
singlet channel with that of the leading-order thermal perturbation theory. Similar to the case of mofiiT), we find 
that the lattice results of mo.2{T) are much larger than the prediction of the thermal perturbation theory at the 
leading-order. 

In Fig. [331 (right), we plot the lattice results for the ratio m£i.2/"^D,o and compare them with ([5^ . We find that 
this ratio also deviates from the prediction of the leading-order thermal perturbation theory. We note that, with 
the p4- improved staggered quark action, the ratio ^2/^-0,0 was reported to agree with S/Stt^ at T ^ 1.5Tpc [63{. 
Similar discrepancy between Wilson and staggered type quark actions has been already reported for Debye screening 



masses at /i, = [2J|. Further investigations at smaller lattice spacings etc. are required to clarify the origin of the 
discrepancy. At /i^ = 0, it was shown that the discrepancy with the thermal perturbation theory is largely removed 
for mD,o{T) with the improved Wilson quark action when we include the next-to-leading-order contributions [23]. 
Thus, a higher order calculation of the thermal perturbation theory at finite will also be important to understand 
the results obtained on the lattice. 
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FIG. 33: (Left) mD,2{T) for the color singlet channel. Dashed lines represent the prediction of the leading-order thermal 
perturbation theory for k = nT, 2nT and 3nT from above. (Right) raD,2/mD,o for the color singlet channel. The dashed line 
at m-D,ilva-D,Q ~ 3/(87r'^) represents the prediction from the leading-order thermal perturbation theory. 



VI. CONCLUSIONS 



Comparison of results obtained by different lattice formulations is important to estimate theoretical uncertainties 
in lattice QCD calculations. Since most lattice QCD simulations at finite temperatures and densities have been 
performed using staggered-type quark actions so far, studies with a different lattice quark action is particularly 
important. In this paper, we carried out the first calculation of the equation of state at nonzero densities with two 
flavors of improved Wilson quarks. Simulations are performed on a 16^ x 4 lattice along the lines of constant physics 
corresponding to mps/mv = 0.65 and 0.80 in the (/?, if ) plane. With Wilson- type quarks, statistical fluctuations of 
physical observables at finite density are much severer than with staggered-type quarks. To tame the problem, we 
combined and developed several improvement techniques. 

Adopting the Taylor expansion method, we calculated the derivatives of pressure with respect to the chemical 
potentials iiq and /i/ up to the fourth order . Using these derivatives, we studied the fluctuations of quark number 
and isospin densities at finite chemical potentials. A quantitative difference between the second derivatives of Xg and 
Xi was observed: Xq shows a peak near Tpc, whose height increases as /i^ increases, whereas xi does not show a clear 
peak near Tpc- These behaviors agree qualitatively with the results obtained using p4-improved staggered fermions, 
and are consistent with the expectation from the effective sigma model. 

With the current statistics, the statistical errors in the results were not small with the simple Taylor expansion 
method. To improve the calculation, we adopted a hybrid method of the Taylor expansion and the reweighting tech- 
niques combined with a Gaussian approximation for the distribution of the complex phase of the quark determinant. 



In a previous study with a staggered-type quark [23[, this method was shown to be efficient to suppress statistical 



fluctuations at finite densities. We found that the statistical errors in the quark number density and the susceptibility 
at finite densities are reduced with the new method. Although the simulations at different temperatures are indepen- 
dent, the resulting T-dependence in the quark number density and the susceptibility turned out to be smooth, and 
the heap in Xq near Tpc became clearer. These results suggest that the sign problem at finite densities is mildened by 
such improvements. 

We also studied the heavy-quark free energies and the Debye screening mass at finite densities in the high temper- 
ature phase. We calculated the Taylor expansion coefficients of the heavy-quark free energies in all color channels up 
to the second order in Hq/T. We found a characteristic difference between QQ and QQ free energies: The inter-quark 
interactions between Q and Q become week, while those between Q and Q become strong, as //^ increases. We also 
calculated the effective running coupling and the Debye screening mass for each color channel up to the second order of 
jjLq. Both quantities show no significant color channel dependence at T ^ 2Tpc. The second order coefficient of Debye 
screening mass, 'mD,2(T), turned out to be positive, implying that the Debye mass becomes larger as increases. We 
note that our mD,2{T) does not agree with the leading-order thermal perturbation theory. Higher orders are required 
to explain the lattice results. 
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Appendix A: Pressure and quark number susceptibility in the free gas limit 



In order to estimate the equation of state at nonzero quark chemical potential, ^ = i^qU in the high temperature 
limit, we calculate the pressure and its derivatives with respect to fj, in the free quark gas limit. Because the effect of 
finite quark mass becomes negligible in the high temperature limit, we discuss only the case of massless quarks. 

The partition function for free Wilson quarks is give by 



Z{K,tj.) = (detM)^f, 



-K 



e^(l - li)5,+i,y + e-^(l + 74)<5,_4,, 



(Al) 



(A2) 



on an x Nt lattice. Note that the clover term vanishes for free quarks. We perform a unitary transformation into 
momentum space (Fourier transformation). 



^fe' = E e-''=^+''^M., ^ UlM^yUyi. 



Here 



-ikx 



1 ^x(i-k) ^ g^^^^ detiU'iU) = det C/^' det C/ = 1 



NfN, 



We then calculate the partition function, 

Z{K,ii) = (detM)^f = (detM)^f, 



-ix{k—l) 



l-i^5]((l-7i)e''*+(l+70e-''0 

-K (e''(l - 74)e''^ + e-''(l + 74)6"''^)]] 
3 

1 — K (2 cos ki — 2i7j sin fcj) — K {2 cos(A;4 — ifx) — 2174 sin(A;4 — i/x)) 



i=l 



where 



27r04 + 1/2) 



ki = 



, j^=0,±l,--- ,Ns/2 for/i=l,2,3 
, j4 = 0,±l,--- ,7Vt/2. 



(A3) 



(A4) 



(A5) 



(A6) 

(A7) 
(A8) 



32 



Introducing a 4 x 4-matrix which is defined by M^i = 6k,iM{k), 

3N{ 



3 

det M(fc)=det I — K (2 cos fc, — 2i7i sin fcj) — K {2 cos(fc4 — ifi) — 2174 sin(A;4 — z/x)) 

^3 \2 3 

1 - 2iC ^ cos fci - 2 cos(fc4 -ifJ-)] + 4:K'^ ^ sin^ ki + AK"^ siia^^ki - ifi) 



1 - 8K + 4K^sm 



2 I 



i=l 
3 



4ii'sin2 



i=l 



^4 — l/X 



+AK'^ ^ sin^ fci + AK'^ sm^{ki - ifi) 



sm 



fc4 — z/x 



2 ^ sin' 



2 / 



i=l 



3 



4hE 



sm 



1 sin^ 



2 I k^ — ifi 



\ i=l 

(det(aoJ + aini + 02*72 + asijs + aii^i) = {oq + + + «3 + "I)^) 
In the massless quark Hmit K = 1/8, 

16 " 



sin^ /ci 



det M(fc) 



+ B\k) + 4(S(fc) + 1) sin^ 



fc4 — j/i 



where 



A{k) = ^sin^ fci, = 2^sin2 . 

i=i i=i ^ ^ 



We calculate the derivatives of pressure with respect to /x at /x = 0, = 1/8 numerically. 

P 



T4 = 



nW^ ^"^^^'^^ - ^ inz(r = = 0) ) , 



1 d"{p/T*) 



A^/-" a"ln2:(T) 



ju=0 



(A9) 



(AlO) 
(All) 

(A12) 
(A13) 



Here, Z{T,ii) and Z{T = 0,/i) are the partition functions calculated on x Nt and iV^ lattices, respectively, The 
derivative of the normalization ln^(r = 0,/z = 0) in /x is, of course, zero. The derivatives of In 2 at /x = are given 
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FIG. 34: The results of p(/i = 0) (circle), C2 (square) and C4 (triangle) normalized by the values of their continuum limit. 



by 



where 



InZ 
93 InZ 



.3«,|.X:Mdc.Ai,«.6«,X:(|| 



fl2 



= 3A^f 



5^ In det M(fc) = 6A^f 



2?o(fc) 
P3(fc) 



Vl{k) 



P2(fc)Pi(fc) , P?(fc) 



Vl{k) 



'Vl{k) 



6iVf ^ 



^lndctM(fc) 

k 

Vi{k) 



Vo{k) 



vl{k) 



vlik) 



12 



V,{k)Vl{k) Vfik) 



vl{k) 



2?4(fc) 



Po = A{k)+B''{k)+A[B{k) 
2?n:odd = -2i [B(fc) + 1] sinfc4, 
2?«:evcn = -2 [-B(fc) + 1] cosA;4. 



-l]sin2(fc4/2), 



(A14) 
(A15) 
(A16) 



(A17) 



(A18) 
(A19) 
(A20) 



The odd derivatives vanish as in the case of interacting quarks. Since c„ = for n > 4 in the continuum limit, 
we calculate p aX — as well as C2 and C4. The numerical results normalized by the values of their continuum 
Stephan-Boltzmann limit are plotted in Fig. [Mj Circle, square and triangle symbols are the results of p(/Lt = 0), C2 
and C4 for each Nt with Ng/Nt = 4, respectively. The results with Ns/Nt = 8 are also shown by the dashed lines. The 
Ng dependence is found to be negligible. However, the results are much larger than unity for small Nt, suggesting 
sizable lattice discretization effects for Nt < 10. 



Appendix B: Derivatives of In 2 in the Gaussian approximation 



We discuss the error from the Gaussian approximation of a complex phase distribution function. We calculate the 
second and forth derivatives of InZ when the Gaussian approximation is applied, and compare with the exact results. 
Denoting the derivative of In det M as 



N{ 



3" In det M(/^) 



(Bl) 



^J,=o 
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the partition function with the Gaussian approximation Eq. (1371) can be expanded in a power series, 



m 



exp 



F- 



(m=o) 



exp 



(exp(F))i. wo{F) dF 



= / exp 



3! 



2 X (3!)2 



X exp 



2 4 



2 4! 8 8 

We then obtain the second and forth derivatives of InZ at ji — Q. The second derivative is 



wo{F) dF. 



dlnZ 



din') 



InZ 



2 9/i2 



«;o(f) =i((2?2) + (P2)) 



This result is, of course, the same as the exact resuh. Next, we calculate the forth derivative. 



92 InZ 



9(/i2)2 



1 a-^inz 



12 9/^4 



3 12 4 

12 4 



wq{F) dF 



wo{F) dFj 

'{Vl)l ^ {VI) ^{V,); 



woiF) dF~-{{Vl) + {V,)y 



Because F = 272(a*V2) + 0(/) and / (• • ■) p 0[F]wo{F)dF = (• • -OlF]) for any function of F: 0[F] 



{Vl)p {V2)pWo{F) dF = {VIV2) + O(m'). 



(B2) 



(B3) 



(B4) 



(B5) 



Moreover, as discussed in Sec lIVBl Vi is given by a sum of the local number density operator (~ 'ipjo'ip{x)) at 
fiq — 0. If the simulation is performed apart from a singular point, we may adopt the Gaussian approximation for the 
distribution of I?i . In such a case, 2?i satisfies 



1 



'F 3 

Substituting this equation, the forth derivative becomes 
94 InZ 



4 (2?i2?3> + {V4) + 3 (Vl) + (Vt) + 6 {VIV2) - 3 {(vj) + {V-i))' 



(B6) 



(B7) 



This is the same as the exact result. In this calculation, we assumed that the distribution function of the total quark 
number, is of Gaussian at /iq = 0. Within this condition, we find that the Gaussian approximation does not affect 
the calculation of the derivatives of InZ up to 0{^'^). Similar discussion is also possible for the higher order terms of 
/i and one can find out the condition in which the Gaussian approximation is valid for each order of /x. 



Appendix C: Results of expansion coefficients for Ofefi and mo 

To evaluate expansion coefficients of acff and mc, we fit the normalized free energies with (j77p ~ (l78l) . Our results 
of the expansion coefficients together with the quality of the fits are summarized in Tables IVHiXI We adopt the fit 
ranges 0.5 < rT < 1.0 for Eq. (|77p and 0.25 < rT < 1.0 for Eq. (|75|) . These fit ranges are chosen by examining 
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TABLE V: Results ai{T) and x^/Ndf for the fit of 1st order coefficients at mps/mv = 0.65 (left) and 0.80 (right). The 2nd 
parentheses in 3* channel of ai(T) expresses the systematic errors due to difference of the fit range. 



mps/rriY = 0.65 





Qi(r) X 


10 






R^6 


3* 


6 


3* 


1.07 


-0.01(2) 


10.86(501)(42) 


0.36 


1.10 


1.18 


1.64(116) 


3.64(95)(54) 


1.11 


0.52 


1.32 


0.92(66) 


3.58(85)(2) 


2.42 


0.47 


1.48 


1.07(36) 


2.14(35)(21) 


1.27 


0.82 


1.67 


1.49(41) 


0.90(16)(18) 


1.04 


0.86 


2.09 


0.62(13) 


0.54(11)(11) 


0.92 


1.68 


2.59 


0.44(9) 


0.56(10)(7) 


1.91 


2.14 


3.22 


0.26(5) 


0.36(5)(4) 


0.48 


1.25 


4.02 


0.33(5) 


0.29(5)(2) 


1.13 


1.16 



mps/mv = 0.80 





Ql 


(T) X 10 






R = 6 


3* 


6 


3* 


1.08 




3.03(51)(40) 




0.71 


1.20 


1.35(67) 


2.37(44)(68) 


0.98 


1.27 


1.35 


1.35(36) 


1.28(26)(22) 


2.18 


0.88 


1.69 


0.65(10) 


0.92(14)(13) 


1.05 


1.21 


2.07 


0.50(7) 


0.36(6) (16) 


1.87 


2.81 


2.51 


0.45(7) 


0.38(4)(1) 


0.70 


0.38 


3.01 


0.23(3) 


0.34(3)(1) 


1.83 


2.04 



TABLE VI: xV^DF for the fit of 2nd order coefficients at mps/my = 0.65 (left) and 0.80 (right). 



mps/rriy = 0.65 


mps/mv = 0.80 


T/Tpc 


R=l 


8 


6 


3* 


T/Tpc 


R=l 


8 


6 


3* 


1.07 


0.87 


1.29 


0.64 


1.05 


1.08 


0.63 


1.15 




1.82 


1.18 


0.43 


0.85 


0.64 


1.04 


1.20 


1.00 


1.70 


2.17 


1.46 


1.32 


1.86 


0.95 


1.17 


2.22 


1.35 


0.99 


0.65 


0.64 


0.46 


1.48 


1.02 


1.56 


1.10 


1.32 


1.69 


2.83 


2.49 


1.53 


1.44 


1.67 


1.12 


1.73 


1.22 


0.61 


2.07 


0.95 


0.64 


1.05 


0.98 


2.09 


1.01 


0.96 


2.43 


1.84 


2.51 


1.83 


1.28 


0.73 


1.25 


2.59 


1.19 


1.49 


1.66 


1.21 


3.01 


1.08 


1.76 


1.08 


0.59 


3.22 


1.02 


1.83 


1.98 


0.90 












4.02 


1.61 


0.72 


1.10 


1.86 













the fit range dependence as follows. Let us denote the fit range as i?ini < rT < Ran- We find that the fit results 
are insensitive to i?fin when Run is sufficiently large. To evaluate the sensitivity on i?ini, we introduce i?ini+2 as the 
next-neighboring longer distance on the lattice. For example, when — 0.5 at Nt = 4, the lattice distance of the 
point is 2 and the next- neighboring longer distance is Vl^ + 1^ + 22 = y^, and thus i?ini+2 = V6/4. Similarly, when 
^ini — 0.25 at Nt — 4, i?ini+2 — Then, we estimate the systematic error due to the fit range by the difference 

of the fit results between and i?ini+2 with fixed i?fin- The systematic errors are shown in the second parentheses 
for ai{T) of color-antitriplet channel in Tab. fVl for m£i_2(T) of color-singlet channel in Tabs. IVllIl and HXl We find 
that the systematic errors are almost comparable with the statistical errors, except very close to Tp^. 
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TABLE VII: Results of 02 (T) x 10 at mps/my = 0.65 (left) and 0.80 (right). 







nips/mv = 


0.65 




mps/my = 0.80 


T/Tpc 


R=l 


8 


6 


3* 




R=l 


8 


6 


3* 


1.07 


0.05(67) 


0.40(365) 


-0.03(5) 


-3.83(346) 


1.08 


-0.09(44) 


-0.08(15) 




-0.06(60) 


1.18 






— • oy yo i-O J 




1.20 


-0.30(34) 


-8.42(1002) 


-1.93(196) 


1.03(42) 


1.32 


0.55(71) 


0.31(91) 


-0.80(161) 


-1.20(94) 


1.35 


0.21(22) 


0.01(71) 


-0.81(58) 


0.18(25) 


1.48 


0.35(34) 


-0.70(216) 


-1.66(151) 


-0.41(59) 


1.69 


-0.08(13) 


0.38(35) 


-0.04(28) 


0.06(18) 


1.67 


0.03(24) 


-1.91(127) 


-1.15(97) 


0.27(28) 


2.07 


0.02(7) 


-0.55(31) 


0.01(22) 


-0.16(9) 


2.09 


-0.16(13) 


-0.05(56) 


-0.13(34) 


-0.24(15) 


2.51 


0.03(6) 


-0.24(20) 


-0.08(15) 


0.01(7) 


2.59 


0.07(11) 


-0.03(33) 


-0.42(22) 


-0.21(15) 


3.01 


0.01(3) 


0.00(10) 


-0.08(7) 


-0.07(5) 


3.22 


-0.04(8) 


-0.20(24) 


-0.18(15) 


0.22(9) 












4.02 


-0.04(7) 


-0.21(22) 


-0.28(13) 


-0.01(7) 













TABLE VIII: Results of mD,2{T) at mps/my = 0.65 (left) and 0.80 (right). The 2nd parentheses in the color-singlet channel 

expresses the systematic errors due to difference of the fit range. 





//)ps/"'\" = 


().") 






nil 


's/lllY = 


.80 




T/Tpc 


R=l 


8 


6 


3* 


T/Tpc 


R=l 


8 


6 


3* 


1.07 


1.47(26) (52) 


0.01(348) 


2.17(229) 


1.03(42) 


1.08 


1.17(15)(10) 


1.03(221) 




1.22(20) 


1.18 


0.55(19)(26)- 


-3.18(180) 


-0.27(71) 


0.58(23) 


1.20 


0.76(13)(12)- 


-0.41(69) 


0.63(31) 


0.99(13) 


1.32 


0.84(18)(28)- 


-0.19(100) 


-0.20(57) 


0.49(19) 


1.35 


0.69(8)(1) 


0.40(37) 


0.26(17) 


0.68(8) 


1.48 


0.73(13)(2) 


0.27(72) 


0.06(48) 


0.55(19) 


1.69 


0.40(7)(31) 


0.73(24) 


0.50(17) 


0.61(9) 


1.67 


0.40(18)(35)- 


-0.31(55) 


0.10(27) 


0.53(19) 
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0.45(7)(12) 


0.30(13) 


0.43(11) 


0.37(7) 
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0.38(8)(22) 
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0.14(14) 


0.33(10) 


2.51 


0.38(6)(5) 


0.27(11) 


0.40(9) 


0.39(7) 


2.59 


0.41(7)(27) 


0.47(18) 


0.13(9) 


0.24(9) 


3.01 


0.37(4) (7) 


0.31(8) 


0.29(6) 


0.35(4) 


3.22 


0.44(8) (3) 


0.30(16) 


0.30(12) 


0.65(10) 












4.02 


().11(1())(2<S) 


().ii(il)) 


0.20(12) 


0.i:!(8) 
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TABLE IX: Results of nioMT) determined with the assumption a2{T) = at mps/mv = 0.65 (left) and 0.80 (right). The 
2nd parentheses in the color-singlet channel expresses the systematic errors due to difference of the fit range. 





mps/mv ~ 


0.65 
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8 


6 


3* 




7?= 1 


8 


6 
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1.33(19)(10) 


0.81(153) 


3.09(166) 


1.43(32) 
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0.45(11) 


0.67(6) 
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0.49(11) 
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0.40(5) 


2.59 


0.39(5)(6) 


0.54(10) 
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